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Abstract 


Using  the  restricted  three  body  problem,  the  equations  of  motion  (EOM)  and 
Hamiltonian  are  computed  for  die  moon's  orbit  in  pt^ical  variables.  A  periodic  orbit  is  found 
in  die  vicinity  of  the  moon's  orbit,  and  classical  Floquet  theory  is  applied  to  the  periodic  orbit  to 
give  stability  informatitMi  and  the  complete  solution  to  the  equations  of  variation.  Floquet 
theory  also  supplies  a  transformation  from  pl^ical  variables  to  modal  variables.  This 
transformati(Ni  to  modal  variables  is  made  canonical  by  constraining  the  initial  transformation 
matrix  to  be  symplecdc.  Actual  lunar  data  is  used  to  calculate  the  modes  for  the  real  moon's 
orbit.  Once  satisfied  diat  die  moon's  real-world  modes  are  in  (or  near)  the  linear  regime  of  the 
periodic  orbit,  the  modal  EOM  are  found  by  doing  a  perturbation  expansion  on  the  new  modal 
Hamiltonian.  The  modal  results  from  die  real  lunar  orbit  are  compared  with  the  modal 
EOM/expansion  results.  The  modal  expansitxi  proves  to  be  an  accurate  solution  to  the  moon's 
orbit  given  enough  expansion  terms. 
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THE  MODAL  SOLUTION  TO  THE  MOON'S  ORBIT  USING 
CANONICAL  FLOQUET  PERTURBATION  THEORY 

I.  Introduction 

The  study  of  die  Lunar  Theory,  which  is  essentially  a  particular  case  of  the  three  bocfy 
problem,  began  in  1687  widi  the  introduction  of  Sir  Isaac  Newton's  Principia.  Significant 
contributions  were  made  over  die  last  three  c^ituries  including  those  by  D*  Alembert,  Euler, 
and  Laplace,  however  most  of  these  approaches  are  implicitly  based  on  two  body,  K^lerian 
orbits  as  reference  trtyectories  for  classical  analytic  perturbation  tediniques.  The  first  real 
departure  fiom  this  classical  Lunar  Theory  wasn't  until  1877  vidien  G.W.  Hill  expanded  on 
Euler's  three  body  problem  by  choosing  rectangular  coordinates  that  rotate  with  the  Sun's  mean 
angular  velocity  instead  of  the  moon's  mean  angular  velocity.  In  addition.  Hill  made  some 
simplifications  to  vdiat  is  now  known  as  the  restricted  three  body  problem.  The  result  of  Hill's 
work  is  a  periodic  reference  orbit  in  the  three  body  problem  coordinates. 

One  hundred  years  later,  in  1981,  Wiesel  demonstrated  die  concept  of  applying  the 
classical  Floquet  problem  to  the  periodic  reference  orbit  to  get  a  Floquet  mode  reference 
solution.  The  advantage  of  diis  over  die  classical  refermce  solutions  is  that  the  Floquet  modes 
already  include  some  perturbing  diaracteristics  (such  as  precession  of  the  orbit  plane)  in  the 
reference  soluticm,  before  the  perturbation  analysis  is  done.  More  recently,  Wiesel  and  Ptdden 
(1992)  improved  on  the  concept  by  constraining  the  Floquet  transformation  to  modal  variables 
to  be  canonical.  The  result  is  a  camxiical  Floquet  referoice  solutitm  that  is  quite  applicable  for 
the  Lunar  Theory.  However,  while  Pdden  (1992)  did  illustrate  its  utility  with  a  Sun- Jupiter 
restricted  diree  body  example,  diis  referoice  solution  tedinique  has  not  yet  been  validated  widi 
real  celestial  data  or  for  a  diree  dimoisional  case. 
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The  purpose  of  this  study  is  to  construct  the  three  dimensional  application  of  the 
restricted  three  body  canonical  Floquet  reference  solution  to  a  celestial  body  for  which  we  have 
a  plethora  of  historical  data;  the  moon.  Specifically,  we  will  find  a  periodic  orbit  in  the  vicinity 
of  the  mo(Mi  (Hill's  orbit),  get  the  Floquet  nKxles,  canonicalize  them,  constrain  the  dynamics  to 
match  die  real  data,  and  finally  get  die  modal  reference  solution.  This  work  accomplishes  three 
things; 

1 .  It  provides  a  second,  dynamically  more  chall^ging  system  to  demonstrate 
canonical  Floquet  theory. 

2.  It  investigates  the  applicability  of  the  canonical  Floquet  theoiy  to  an  actual 
system  that  is  constrained  to  obey  historical  data. 

3.  It  lays  the  foundation  for  applying  further  analytic  perturbation  techniques  to 
the  moon's  modal  solution  as  a  reference  tr^ectoiy. 

Chapter  Two  identifies  the  specific  background  work  that  was  done  prior  to  this  study,  and 
explains  the  differoices  between  diis  and  previous  applications  of  Floquet  reference  solutions. 
Chapter  Three  is  a  detailed  discussion  of  the  theory  used  to  investigate  the  modal  solution, 
vsdiile  the  actual  software  used  to  implement  die  dieoiy  is  described  in  Chapter  Four.  Chapter 
Five  is  the  results.  It  includes  an  analysis  of  die  modal  solution  as  well  as  real  moon's  linearity 
with  respect  to  the  periodic  orbit. 
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11.  Historical  Development 


2.1  Dynamics 

The  groundwork  for  the  restricted  diree  body  problem  was  laid  out  by  Euler  in  1772. 

In  applying  this  restricted  three  body  problem  to  the  Lunar  Theory,  Hill  made  some 
simplifications  that  give  us  the  foundation  for  our  dynamic  model.  Plenty  of  information  was 
available  on  this  subject  (Brouwer  and  Clemoice,  196L336  or  Szebehely,  1%7:602). 

2.2  Floquet  Reference  Solutions 

Since  1883,  classical  Floquet  analysis  of  periodic  systems  was  a  common  technique 
only  in  so  far  as  determining  stability  informaticxL  The  idea  of  expanding  the  utility  of  Floquet 
theory  to  get  perturbation  reference  solutions  for  periodic  orbits  was  introduced  by  Dr.  Wiesel 
in  1980.  As  mentioned  in  Chapter  One,  Dr.  Wiesel  that  went  on  to  develop  Floquet  nnode 
reference  solutions  specifically  for  the  Lunar  Theory  in  1981.  This  Floquet  perturbation  theory 
for  periodic  orbits  was  fiirdier  demonstrated  by  Ross  (1991).  His  work  was  actually  the  first  to 
apply  die  classical  Floquet  perturbation  theory  to  a  system  using  the  restricted  three  body 
problem  for  the  dynamics.  The  system  Ross  chose  was  the  Sun-Jupiter  tystem. 

2.2  Canonical  Floquet  Theory 

The  classical  Floquet  problem  naturally  introduced  a  transfimnation  from  idiysical 
variables  to  modal  variables.  The  next  st^  in  the  development  of  this  theory  was  made  by 
Wiesel  and  Pdilen  (1992)  in  realizing  diat  this  Floquet  transfixination  needed  to  be  canonical. 
In  dieir  paper,  they  presented  die  mediod  for  making  the  transformation  canonical  which 
enabled  canonical  perturbadtm  theory  to  be  applied  properly  to  Hai:!iiltonian  systems.  Finally, 
Pohloi  (1992)  illustrated  this  canonical  Floquet  perturbation  dieory  Ity  applying  it  to  the  same 
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restricted  three  body  problem  done  Ross  (1991).  Like  Ross,  Pohlen  used  the  Sun-Ji^iter 
system  because  its  high  mass  ratio  made  it  a  good  diagnostic  case. 

2.4  Application 

The  restricted  three  body  system  used  by  Ross  and  Pohlen  was  confined  to  two 
dimoisions.  Because  of  the  moon's  out-of-plane  motion,  the  problem  in  diis  work  is  ejqMUided 
to  three  dimensions.  Also,  in  Pohlen's  study,  an  effort  was  made  to  transform  complex 
eigensystems  into  real  eigensystems  to  yield  real  modal  vectors.  In  this  study,  we  allow  die 
eigoisystems  and  resulting  modal  vectors  to  be  complex.  Finally,  i^hore  Pohlen  used  die 
restricted  three  body  Sun-Jiqiiter  model  as  die  exact  representation  of  the  modal  solution,  we 
go  on  to  constrain  the  Sun-Eardi-Moon  model  to  obey  actual  moon  ephemerides. 
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Ill  Theory 


3.1  Dynamics 

The  first  stq)  to  finding  the  complete  modal  solution  is  to  decide  on  die  c^namic  modd 
to  be  used  for  d^  systeni  In  this  case,  the  system  consists  of  the  eardi  and  moon,  the  sun,  and 
other  lower  order  perturbing  contributors.  The  restricted  three-body  problem  (R3B)  is  used  to 
define  the  dynamics  of  the  earth-sun-moon  system  Perturbation  theoiy  will  account  for  the 
omitted  terms. 

3.  J.  I  The  Restricted  Three-Botfy  Problem 

The  general  restricted  three-body  problem  was  first  introduced  by  Leonard  Euler  in  his 
memoir  on  his  second  lunar  dieoiy  in  1772.  The  definition,  according  to  Szebehdy,  is  this: 

Two  point  masses  m,  and  called  the  primaries  revolve  around  their 
center  of  nuiss  in  circular  orbits,  hi  the  plane  of  their  motion  moves  a 
third  body  with  infinitesimal  mass,  not  influencing  the  motion  of  the 
primaries.  Assuming  Newtonian  gravitational  forces,  find  the  behavior 
of  the  diird  borfy.  (Szebdiely,  1967:557) 

In  diis  study,  we  assume  the  conditions  above  are  ai^licable  where  the  two  primaries 
are  die  earth  and  sun,  and  the  body  of  interest  (the  third  bocty)  is  the  mooa  We  assume  that  the 
earth  and  sun  are  both  point  masses  and  significantly  more  massive  than  the  moon,  and  diat  die 
moon  has  a  negligible  afiect  on  die  earth  and  sun.  What  we  now  have  is  the  basic  setup  fi>r  G. 
W.  Hill's  approach  to  lunar  dieory.  Hill  tordc  die  geimral  restricted  three-bo^  problem,  ^plied 
the  earth-sun-mooi  system,  and  made  diree  sinqilificatioos  to  the  system; 

1.  The  solar  parallax  is  zero.  Hie  disturbing  function  is  truncated. 

2.  The  lurun  inclination  is  zero. 

3.  The  solar  eccentricity  is  zero. 
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Further  details  of  Hill's  method  am  be  found  in  either  Brouwer  and  Clonence  (1%1  ;336)  or 
Szebehely  (1967:602). 

Before  foe  a)ordmate  syston  for  this  problon  is  established,  we  need  to  introduce 
some  dimensionless  variables  mass,  length  and  time: 


m, 


A/, +A/j  '  A/,+Afj  '  A/, 

5, 

j,  = - = —  j,  =  — = — 

‘  *  5, +5j 


(1) 


where  A/, ,  A/j .  and  A/j  are  foe  masses  of  the  sun,  earth,  and  moon  respectively.  and 
are  foe  distances  of  foe  sun  and  earth  from  the  center  of  mass  of  foe  system.  The  total 
dimensionless  distance  between  foe  earth  and  sun  is  set  equal  to  one,  as  is  the  total 
dimensionless  mass  of  foe  two  primaries. 

j, +^2=1  m, +OTj=l  (2) 

Using  foe  definition  for  center  of  mass  position  ( 5,  )  wifo  respect  to  m, , 

Em, 5,  m,  .0 +m2(s,  +  )  _ 

Erw,  m,  +mj 

we  can  now  define  all  four  dimensionless  variables  in  terms  of  one  non-dimensional  paramefor. 


5,  =m2  sfi 

Sj=w,=l-^  (4) 

The  orbital  period,  r, ,  of  foe  earth  around  foe  sun  is  also  defined  using  the  dimensionless 
parameters  as 


■  VG(m,+m2 


) 


(5) 


where  G  is  foe  universal  gravitational  constant  If  G  is  set  to  one,  t^  =  In  .  In  addition,  the 
angular  velocity,  m ,  of  the  primaries  about  the  center  of  mass  is 
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A  synodic  coordinate  system  is  used  as  deBned  by  Szebehely  (1967:9)  where  the  origin 
is  the  center  of  mass  of  the  two  primaries.  The  X>axis  is  a  line  between  the  sim  and  eardt  The 
sun  is  a  distance  n  from  the  origin  in  the  positive  X  direction,  while  the  earth  is  ( l-/i )  in  the 
negative  X  directioa  The  Y-axis  is  such  diat  the  XY  plane  defines  the  plane  in  which  the 
primaries  rotate.  The  coordinate  system  is  synodic  since  the  X  and  Y  axes  rotate  mih  the 
primaries  around  die  center  of  mass/origin  (see  Figure  1).  The  vectors  r,  i;,and  Fj  are  the 
position  vectors  of  the  moon  widi  respect  to  the  oentet  oi  mass,  the  sun,  and  the  earth 
respectively. 


Figure  1.  The  Restricted  Three  Body  Coordinate  Syst«a 
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This  coordinate  syst«n  is  the  same  as  described  Ross,  with  die  addition  of  the  Z- 
axis.  The  addition  of  a  third  dimeiBion  has  also  been  discussed  in  Szebehely  (1967;SS7). 

3. 1. 2  The  Restricted  Three-Btxfy  Equations  of  Motion 

Let  the  generalized  coordinates  (q^)  «id  momenta  )  be  defined  as 

qt=x.y,z  Pi=Ps.Py.P. 

Then,  the  position  of  the  moon  with  respect  to  the  center  of  mass  is 


r  = 


\^J 


A  tK 


or  r  =  xX+yY  +  zZ 

and  die  inertial  velocity  of  the  moon  in  terms  of  X ,  Y ,  Z  coordinates  is 

V  =  iF  +  (ij?xr)  = 


'x^ 

y 

+ 

X 

Jj 

(7) 


^x-y^ 
y  +  x 

\  i  / 

Using  Eqs  (7)  and  (8),  the  kinetic  (T)  and  potential  (  V)  energies  of  the  moon  are 
calculated  as  well  as  the  specific  Lagrangian  (/,): 

1  _  _  1 


(«) 


7  =  .  V  =  Im, [(i  -yf+(y  +  xY+z^] 


(9) 


Jf  _  -Gniimj  ^  -Gm^ntj 

Kl 


-^3(1-//)^  -WjU) 


(10) 
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I*  sT-V  =  +(j>+*)*  +j^j  +  --  +  . 

/iij  2  ‘  •'  'i 


where 


rx=M  =  ^{x-Mf  +y^ 


r,=\f,\  =  ^ix  +  l-MT+y^-i-z^  (12) 

Taking  appropriate  partial  differentials  of  the  Lagrangian  gives  expressions  fix’  the  generalized 
momenta  in  terms  of  their  conjugate  generalized  velocities  and  vice  versa: 


dl  . 

x  =  p,+y 

(13) 

dl  . 

y  =  P,-x 

(14) 

dl  . 
dz 

z  =  p. 

(15) 

The  Hamiltonian,  X,  is  now  defined  as 

X  =  (« 

only  after  eliminating  die  ^/s  infiworofdie  /r/s.  Substituting  Equations  (11)  and  (13-15) 
into  Equation  (16)  yields 

^  =  ^Py  +  (!' 


r,  r. 


To  get  the  equations  of  motion  for  die  state  vector 

=[x.y,z,p,.p,,p.] 


9 


we  lee  the  Hamilton  canonical  equations: 


.  _dX  .  dX  .  _ ,  ,  , 

a-  ’  t  -  1.2,3 
oPi  oq, 

Thwefore,  the  equations  of  motion  for  this  problon  are 

.  dX 

.  eH 

y  =  ^=Py-X 

^Py 

.  dX 


p{x  +  l-p) 


U-M)'y  py 


dX  U'Z 


-pyz  p. 
r*  ~ 

M  '2 


where  the  p/s  are  the  inertial  veiodty  components  resolved  on  the  rotating  frame. 
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3.2  Periodic  Orbits 

In  the  last  section  we  established  the  equations  of  motion  (EOM)  for  the  restricted 
thi  ee-body  problem.  Now  we  need  to  detomine  the  initial  conditions  that  create  a  periodic 
(Mbit  (A^ch  generally  describes  the  moon's  motitm.  Obviously,  the  moon  is  nor  in  a  perkxlic 
orbit  for  the  resticted  problem's  cocMdinate  system  Because  of  eccentricity,  inclination,  the 
precession  of  the  plane  of  the  (Mbit,  and  many  other  perturbations,  the  lunar  orbit  does  not  (XMne 
back  to  meet  itself  every  month ...  but  it  comes  dose.  The  supposition  is  that  it  comes  close 
enough  to  at  least  be  in  ((m  near)  the  linear  region  of  a  periodic  (Mbit 

A  perkxlic  (Mbit  is  (xie  which  closes  on  itsdf  after  each  revolution.  In  (Mher  wtmcIs,  tlto 
initial  state  vector  (of  a  body  in  a  perkxlic  (Mbit)  will  be  equal  to  the  state  vector  at  any  multiple 
of  the  period: 

x(0)  =  Jc(r)  (20) 

However,  given  uncertainties  in  the  clynamical  nxxld,  it  is  difficult  to  determine  what  initial 
conditirNis  should  be  used  to  create  the  peri(xlicity.  In  practice,  once  a  set  of  initial  ccxiditions 
has  been  chosen  and  die  orbit  integrated,  one  wiU  firxl  that  the  initial  and  final  (xxiditions  will 
not  agree."  (Ross,  1991 : 10).  It  is  useful,  dierefixe,  to  integrate  the  equatkms  of  variatkxi 
(EOV)  al(Hig  widi  the  EOM,  because  the  equati(xis  of  variatirxi  allow  us  to  handle  nearby 
orbits.  It  is  due  to  these  equations  of  variation  that  we  can  correlate  differences  in  the  boundary 
(XMiditirxis  to  correctirxis  in  the  initial  conditions.  As  will  be  discussed  shortly,  diis  sets  up  an 
obvious  iteradrxi  for  crxiverging  on  the  correct  initial  conditions  diat  satisfy  Eq  (20)  and  create 
die  desired  perirxlic  orbit 
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3.2.  J  Equations  of  Variation 


In  order  to  iterate  on  the  proper  initial  conditk»s  required  for  a  p^iodic  orbit,  it  is 
necessary  to  integrate  the  equations  of  variation  as  well  as  the  equations  of  motion.  The  EOM 
can  be  written  as 

x  =  Z^  =  /{x.t)  (21) 

OX 

where  Z  is  die  correlation  matrix  which  has  die  form 

Z  =  j  j  I  =  identity  matrix  (22) 


This  correlation  matrix  follows  the  identities:  =  Z  '  =  -Z  .  If  we  define  the  state  ( x  )  as 

the  periodic  trajectory  ( )  plus  a  small  variation  (  ^x  ), 


x  =  Xp+<yx  (23) 

substitute  into  Eq  (21)  and  eiqiand  in  a  Taylor's  series  about  ^x  =  0  (or  x  =  Xp  ),  we  get  die 
equations  of  variation  (Wiesel,  1993:1 14) 


=  A{t)dx  (24) 

where  ^4  is  a  square  matrix  of  partial  derivatives  of  the  equations  of  motion  with  respect  to  the 
state  variables,  evaluated  on  the  periodic  trtgectory.  It  is  a  fiinction  of  time  only.  For  this  sturfy 
diei4-matrix  is 


12 


dx  dx  dx 
dy  dz  dp^ 


dx  dx 

^P,  dp. 


A{f)  = 


0  1  0  10  0] 


-1  0  0  0  1  0 


0  0  0  0  0  1 


Am  >!«  >4,3  0  1  0 


■^51  Af2  A^i  10  0 


^  1^.  ^  ^  0  0  OJ 


where 


dp.  .  i^-p)  .  M  .  3/i(x  +  l-Ai)^ 

=  — —  =  1 - z —  + - - - r  + - i - 


dx 


ri 


A  ^Py  _  ,  p  .  3ai/ 


r,  r, 


dp,  _  3(l-/i)z*  p  Zpz 

,  T  m  ""  «  ^  C 


ar 


3  5 

'2  r: 


'^SI  “  A,i2 


ap,  3(l-/<)(x -/<);;  3;3(x  +  1-/i);' 

_  c  '■ 


•^<1  “  A^j 


_  ap,  ^  3(l - /<)(x - p)z  ^  3p(x  +  l-p)z 


dPy  3{\-p)yz  Zpyz 


-^62  "  ^J3  “  a.  “  J  ^  5 

ar  r,  r. 


We  now  introduce  die  square  matrix,  0(/,/o) ,  which  is  made  up  of  N  =  6  columns 
(each  of  wMdi  independendy  satisfies  Eq  (24)).  <li(^  O  >s  called  die  smi^  transitkm  matrix 
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and  it  maps  changes  in  die  initial  conditions  to  changes  in  the  final  cotiditions  through  the 
relationship 

5Jf(/)  =  <!>(/./„  (27) 

4>  therefiire  satisfies  the  equations  of  variation  which  now  take  the  fimn 

<i>(/./o)  =  i4(/)0(/./J  (28) 

where  die  initial  condition,  ) ,  equals  dtt  identity  matrix  ( I ). 

3.2.2  Our  Periodic  Orbit 

We  now  have  the  necessaiy  tools  to  find  die  periodic  orbit  fix  our  system.  The  s^up 
for  our  periodic  orbit  is  shown  on  Figure  2.  It  is  a  symmetric  periodic  mbit  with  zero 
inclination,  therefore  we  set: 


x{o)  =  specified  />,(o)  =  0 

y{0)  =  0  Pj^{0)  =  ^cified  (29) 

z(0)  =  0  a(0)  =  0 


Aldiougb  the  moon's  motkm  is  not,  in  reality,  cmnpletely  confined  to  the  plane  of  the  primaries' 
motion,  we’ll  let  it  be  so  fiir  the  periodic  orbit  and  use  other  techniques  to  handle  the  out-of- 
plane  motioa  This  assunqition  was  also  made  by  Hill. 

The  period  of  die  moon  in  the  restricted  problem  is  calculated  firom  a  real  ephemeris  as 


x=2n- 


=  2jt 


29.530589 


(30) 


365.256363 

vdiere  r  is  the  paiod  of  the  sidereal  year  and  is  the  mean  period  firom  new  moon  to 


new  moon  as  tabulated  in  l3aieAstroneaiti<xil Almanac  (U.S.  Naval  Obsovatory,  1986:C1T32). 
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Figure  2.  Periodic  Orbit  in  the  Restricted  Three  Bo(fy  Problem. 


At  this  point,  we  also  determine  a  value  for  the  parameter  /r : 

M 

^  =  —  =  3. 00348069  x  10^  (31) 

^  '  A/, '  ’ 

Here,  A/,  =  1  ,  and  A/j  =  1/332,946.038  solar  vaasses  (U.S.  Naval  Observatory,  1986:K6). 

The  general  ^proach  is  to  start  with  a  guess  for  die  initial  conditions  based  on  two> 
bo(fy  dynamics  and  iteratively  integrate  the  state  vector  and  matnx  forward  in  time  to  one 
period  using  the  EOM/EOV  until  Eq  (20)  is  satisfied..  Specifically,  we  want  y(0)  =y{T)  =  0 
and  Px(0)=p^(i)  =  0  .  zandp^^sU^atzeroby  thonselves,  and  z  and  p^  will  return  to 
their  initial  values  due  to  symmetry  as  long  as  y  and  return  to  zero.  After  each  int^ratkm, 
the  error  at  t  =  r  is  simply; 


IS 


So  die  correction  to  die  final  state  is  just  -e  .  Now  we  want  to  know  how  to  change  the  initial 
conditions  to  effect  die  corrections  in  the  final  conditions.  Using  Eq  (27),  and  removing  all  of 
the  superfluous  elements,  we  have; 

Sx{  r)  =  <1>(  r,0)^jc(o) 

^25 

By  solving  for  ^Jc(o) ,  we  find  the  correctitxis  to  the  initial  conditions.  Those  corrections  are 
added  and  the  process  is  started  again.  This  is  done  iteratively  until  some  tolerance  is  met  (we 
used  10'°). 


^x(O) 

SPy{0) 


(33) 


3.2.3  First  Guess  for  Initial  Conditions 

To  get  our  first  pass  guess  for  initial  eruditions,  we  used  two  body  dynamics  solutions. 
As  shown  in  Figure  3,  the  initial  x  position  in  the  restricted  diree  body  coordinates  conqionent 
is 


=- 


(l-p)- 


/HOOtt  I 

lAU  J 


Xo=- 


(1-p) 


385,000  km 


1.4959787  x  10*  km^ 
wiiile  the  initial  y  velocity  component  (in  R3B  coordinates)  is 


=  -.997423  AU 


(34) 


yo  = 


I  Pe 

^^moon 


3.986012  X 10^ 

- 2^- =  1.0175  — =.034162  — 

°  385,000  km  s  TU 


(35) 
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v/here  is  the  gravitational  paramo  the  earth  ukI  is  the  approximitte  niMn 
radius  of  the  moon  with  respect  to  earth's  center.  But  since  we  need  the  momentum  in  the^ 
direction,  we  use  the  equations  of  motion  to  get 

Px,  =  A =  --963261  (36) 

This  gives  us  die  first  guess  for  the  initial  conditions  which  looks  like  this; 

x(0)  =  [-.997423  0  0  0  -.963261  of  (37) 

After  doing  the  iterative  correcdons,  the  final  result  is  the  initial  condition  state  vector  for  our 
periodic  orbit 

jCp(0)  =  [-  997456  0  0  0  -.965393  O]^  (38) 


Figure  3.  Two  Body  Initial  Guess  in  die  Restricted  Three  Body  Coordinates. 
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3.3  Classical  Floquet  Theory 

The  basic  contribution  of  classical  Floquet  theory  is  that  the  periodic  state  transition 
matrix,  0 ,  can  be  writtoi  as 

=  (39) 

where  F  is  a  periodic  matrix,  and  7  is  a  constant  matrix  of  system  frequencies  in  Jordan 
normal  form.  These  system  frequencies  are  called  Poincare  oqxxients.  If  we  can  find  the 
constant  J  matrix,  as  well  as  the  F  matrix  over  one  period,  then  we  have  <I>  over  oiw  period 
and,  therefore,  we  have  <I>  for  all  time. 

To  determine  J  and  F ,  the  equations  of  variation  (Eq  (28))  are  integrated  forward  in 
time  to  one  period  (  r).  Conveniently,  this  was  already  done  in  the  last  iteration  of  finding  our 
periodic  orbit  (§3.2.2).  Die  result  is  <I>(  r,0)  ,  called  the  momdromy  matrix,  which  looks  like 

O(r.0)  =  F(ry'F-'(0)  (40) 

But  F(  r)  =  F(o)  because  F  is  periodic,  so 

<p(r,0)  =  F(0)e'^F-'(0)  (41) 

After  rearranging,  we  get 

F'  (0)O(  r.O)F(o)  =  e-'"  (42) 

^^ch  shows  that  F(o)  =  F(t)  is  the  matrix  of  eigenvectors  of  ^(r,0),and  is  the 
diagonal  matrix  of  eigenvalues  of  0(  r.o) .  Therefore,  after  finding  the  mcMKxlromy  matrix, 

die  next  step  is  to  find  its  eigenvalues  and  eigoivectors. 

Given  the  constant  eigenvalues.  A,,,  for  die  system,  thoi  finm  Eq  (42) 

A,  =  (43) 

v\^ere  are  die  system  fi^uencies  (Poincar^  exponents).  Solving  for  (o^  yields 

(0^=—lnX^  (44) 

T 
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which  are  the  constant  diagonal  values  of  the  J  matrix.  These  Poincare  exponents  alwi^ 
occur  as  positive/negative  pairs  for  canonical  systems  (each  pair  corre^xxiding  to  a  mode), 
with  me  exception:  when  =  1  ,  then  there  is  a  repeated  value  of  =  0  .  We  call  this 
pair  of  zeroes  a  degenerate  mode.  For  each  of  these  d^enerate  mode  pairs,  there  is  an  exact 
integral  of  motion,  and  a  value  of  one  appears  in  die  off-diagonal  r^kxi  of  J .  It  is  also 
notewordiy  diat  the  Poincare  exponents  give  important  linear  stability  information  about  dw 
periodic  orbit.  They  "can  be  interpreted  just  like  die  eigenvalues  of  a  constant  coefficioit 
system.  The  imaginary  part  of  o),  is  the  oscillatory  frequency  of  the  mode  r,  while  a  positive 
real  part  indicates  instability"  (Wiesel,  1993: 12S).  For  our  specific  periodic  orbit,  we  get  one 
degenerate  mode  and  two  purely  imaginary  modes.  Because  die  Poincar6  exponents  have  no 
real  parts 

<0,  =  0  +  Of 

ruj  =0  +  0.8853941825307/ 
fij,  =0  +  1.053464567610/ 

0)^  =0  +  0/ 

ojj  =0  -  0.8853941825307/ 

=  0  - 1.053464567610/ 

we  know  the  system  is  stable.  The  monodromy  matrix  <I>(  r,0) ,  initial  eigenvector  matrix 
F(o)  ,  and  J  matrix  are  shown  in  Appmdix  A. 

The  next  step  is  to  get  the  eigmvector  matrix,  F{t) ,  over  (»e  period.  Using  Floquef  s 
Theorem  (Eq  (39))  and  differentiating,  we  find 

So,  by  substituting  into  Eq  (28) 

<i>(/,/J  =  4/)<I>(/,/J 
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fi{l)  +  F(l)j=A{l)F{l) 

Therefore, 

F=AF-FJ  (45) 

is  a  dififerential  equation  for  F{t) .  With  the  initial  condition,  F{0) ,  already  known,  we  can 
int^rate  from  /  =  0  ->  r  to  get  F{t)  over  one  pmod,  which  can  be  reduced  into  a  series  oi 
Fouritf  coefficients  as  described  Brouwer  aKl  Clemence  (1%1 : 108-1 13). 

By  knowing  F{t)  over  one  pafod,  we  also  have  a  coordinate  transformation  from 
pl^ical  variables,  Sx  ,  to  modal  varuAles,  y 

(44) 

which  is  derived  from  equation  (27)  by 

5je(r)  =  F(0e'<-^>F-'0JSx(O 

F-'(l)FX(l)  =  e'<'-^>F-'(l,Mx(lo) 

=  (47) 

Equation  (47)  is  the  modal  coordinate  fimn  of  the  FloquetSolutioa  It  is  a  sdution  to 
the  linear,  constant  coefficient,  modal  equation; 

f  =  Jy  (48) 

w^ete,  for  canonical  systems,  7  has  the  fi)rm  J  =  ZS  ,  S  being  asymmetric  matrix. 

Classical  h  -oquet  theory  gives  us  a  conq)lete  solution  (Eq  (39))  to  the  EOV,  stability 
information,  and  a  transformation  to  modal  coordinates  (Eq  (46)).  To  do  die  perturbation 
expansion  discussed  in  a  later  section,  it  is  necessary  that  the  transfixmation  be  canonical. 
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3.4  Canonical  Floquet  Theory 

In  the  last  section,  we  defined  the  transfisrmation  to  modal  coordinates  as 

We  would  like  diis  modal  transformation  to  be  canonical  so  that  a  new  Hamiltonian  can  be 
derived  fix  the  perturbation  expansion  analysis  in  a  later  sectioa 

The  transformation  above  is  canonical  only  if  F{t)  is  a  symplectic  matrix.  That  is 


Z  =  F^ZF  (49) 

must  hold  for  all  time  v^iere  Z  is  the  correlation  matrix.  This  ensures  that  foe  new  modal 
Hamiltonian,  X{y) ,  follows  Hamilton's  equations  (Eq  (21)) 


j  =  Z 


dX{y) 

dy 


(50) 


where  y  is  foe  modb/ state  vector. 

An  algorifom  was  presented  by  Siegel  and  Moser  (1971 :97>103)  that  nomudizes  F{0) 
to  a  symplectic  matrix,  but  foisalgoiifom  only  accounted  for  a  constant  F  matrix  and  a  non¬ 
degenerate  case.  Later,  Wiesel  and  Pohlen  (1992:6-12)  ejqjanded  the  utility  of  this  algorithm 
to  include  periodic  F  matrices  and  degenerate  modes. 


3. 4. 1  Symplectic  Periodicity  for  Non-Degenerate  Modes 

Because  foe  Poincar6  e^qxments,  alw^  occur  in  positive/negative  pairs  fix' non- 
degoierate  canonical  systons,  the  J  matrix  can  be  written  in  the  fixm 

fn  0  1 

-a] 

vsfoere  is  a  diagonal  matrix  with  one  of  eadi  pair  on  the  diagcxial.  By  dififerentiating 

Eq  (49)  and  recognizing  that  Z  is  constant  matrix,  we  get 

Z  =  F^ZF  +  F^ZF  =  0  (52) 
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Substituting  Eq  (45)  and  Eq  (49)  and  some  correlation  matrix  idenbtm,  this  can  be  rearranged 
to  obtain 

-J^Z-ZJ^O  (53) 

which  proves  to  be  true  since 

^  f  0  -Ql  f  0  nl 

-•'  ^-""=1-0  0  Mo  oh® 

Therefore,  if  F{Q)  can  be  made  symplectic.  then  F{t)  wiU  stay  symplectic  for  all  time  (Wiesel 
andPohlen,  1992.7). 


3. 4. 2  Symplectic  Periodicity  for  Degenerate  Modes 

For  systems  with  a  d^enerate  nnode,  the  J  matrix  will  have  a  pair  zeros  and  take  the 

form 


J  = 


(55) 


where  T*  has  a  1  on  its  diagonal  location  dut  occiq>ies  the  same  row/column  as  the  d^enerate 
mode  and  zeros  eveiy^^me  else.  Just  as  with  the  non-d^enerate  case,  we  substitute  J  into 
Eq  (53)  to  show  dut 


-J^Z-ZJ 


I  0 

0 

si 

l-Q 

-•*r 

(56) 


So,  once  again,  du  same  savts  5)wi^clic,  sarysfymr/Tfeclic  theory  holds  true.  The  next  step  is 
to  actually  make  die  inhid  eigenvector  matrix  s^iqilecdc:  the  subject  discussed  in  the  next 
sectioiL 


22 


3.5  SymplecUc  Normalization 

The  theory  fix'  symplectic  normalization  is  presented  in  detail  by  Wiesel  and  Pohlen 
(1992:8-12),  and  discusses  all  possible  cases.  This  section  only  summarizes  tlw  method  as  it 
applies  to  the  lunar  modes;  that  is,  one  d^enerate  mode  wd  two  non-degen«ate  (pure 
imaginaiy)  modes. 

Right  now,  we  have  the  eigenvector  matrix  {F)of  <t^  \^ch  follows  from  Floquefs 
Theorem.  The  intent  is  toget  another  dgenvector  matrix  (£)whicfa  solves  this  equation 

E-'^E  =  e^^  (57) 

and  at  the  same  time  is  symplectic. 

Z  =  E^ZE  (58) 

We  can  write  the  relationship  of  die  two  eigenvector  matrices  as 

£  =  FD  (59) 

where 


Matrices  </,  and  are  diagonal  matrices  of  multiplicative  scale  factors.  These  scale  factors 
are  arbitrary  with  the  constraint  that  the  rdi  entries  in  d^  and  </,  corresponding  to  repeated 
eigenvalues  (degenaate  mode)  be  the  same.  The  matrix  c  is  zero  everywhere  unless  there  is  a 
degenerate  mode,  in  which  case  there  is  a  1  in  die  rth  entry  along  the  diagonal  corresponding  to 
die  degenerate  mode.  In  our  particular  case,  there  is  a  d^enaate  mode,  so  all  nf  the  caveats 
apply 

Substituting  Eq  (59)  into  Eq  (58)  and  rearranging,  we  get 

D'^ZD'  =  F^ZF  (61) 

where  is  the  inverse  of  the  transpose  of  D .  By  direct  calculation,  die  above  equation 

becomes 
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D'^ZD'  =  ®  =  f^zF  (tt) 

l-d-'d,"  0  J 

Since  we  know  F,  we  can  detennme  the  pioduct  d^dJ  By  settiiig  </,  =  i/j  suggested  by 
Wiesd  and  Pohlen,  we  can  solve  fix  the  diagonal  entries  of  d^  and  d^  by 


With  the  diagonal  values  dPthe  D  noatiix  calculated,  equation  (59)  returns  die 
symplectically  normalized  initial  eigenvector  matrix  E . 

3.5.1  Generalized  Eigenvectan 

As  mentioned  befixe,  our  periodic  orbit  has  a  degenerate  mode,  which  means  one  oi 
the  eigenvectors  (  /^  )  in  the  F  matrix  is  not  linearly  independent  To  continue,  we  need  a 
genoalized  eigenvector  ( ^' )  before  normaiiziog  to  the  F  matrix.  To  get  our  generalized 
eigenvector,  we  need  to  look  at  the  fixm  of  die  coigugate  eigtiiiveclor  fix  the  d^enerate  mode 

/.  =[0  /,  0  /«  0  0]"  (64) 

We  want  linear  independence,  so  if  we  choose  //  to  have  the  fixm 

/;  =  [/,'  0  0  0  /,'  of  (65) 

dien  perpendicularity  between  the  two  is  automatic. 

/.•7;  =  0  (66) 

The  generalized  eigenvector  is  now  calculable  fixxn 

(0-l)f:=Tf,  (67) 

udiwe  r  is  the  period  ofthe  periodic  orbit  and  specific  to  canonical  problems.  Widioutthe 
factor  of  r  included,  the  F  (or  F )  matrix  is  not  periodic. 
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3.6  Modal  Variable  Development 

We  have  tfie  initial  state  vector  for  the  periodic  orbit  in  physical  variables,  jCp(O) ,  as 
well  as  the  initial  symplectically  nmnalized  eigenvector  matrix  £(0) .  Using  equations  o£ 
motion  found  earlier 


X  =Z 

(21) 

E^AE-EJ 

(45) 

the  two  are  integrated  forward  in  time  to  one  period.  As  both  are  periodic,  fo^  can  be  reduced 
to  Fourier  series  coefficients  by  harmonic  analysis.  The  method  for  harmonic  anafysis  is  ^wn 
in  Brouwer  and  Clemence  (1961 :108*1 13),  although  we  use  an  exponential  form  instead  of  the 

sin/cosine  form  to  accommodate  the  complex  valued  E  matrix.  With  these  Fourier 
coefficients,  we  know  Xp(r)  and  E(t)  for  all  time. 

The  next  step  is  to  randomly  select  an  initial  condition  state  vectM-,  Jr(0) ,  for  a  near- 
periodic  orbit  tnyectory,  and  int^rate  forward  in  evenly  ^>aced  time  steps  through  one  period. 
At  each  time  step,  foe  periodic  state  vector  and  E  matrix  are  drived  from  their  Fourier 
coefficients  as 

\{t)=  Zigx)je^ 

j=-m 

E{t)  =  I* 

/=-• 

where  gx  and  ge  are  the  Fourier  coeffidents  and  r  is  the  imagiiuuy  value  The 
variation  in  physical  variables  Sx{t)  is  also  calculated  at  each  time 

«(/)=*(/)-*,(»)  («() 

vfoich  is  finally  transfi>rmed  to  modal  variables  using  the  now  fomiliar  rdation 

y{t)  =  E-'{t)Sx{t)  (69) 

Rem^nber  that  the  E  matrix  replaces  the  F  matrix  fi>r  a  canonical  transfiirmation. 
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We  now  have  the  modal  coordinate  fixm  of  the  Floquet  solution  to  any  orbit  that  is  near 
our  periodic  orbit.  The  firnn  of  these  modal  variables  is  special.  For  the  non-degenerate 
modes,  an  unususal  symmetiy  exists  between  the  conjugate  modal  variables,  and  ; 
their  real/imaginaiy  parts  are  switched  and  of  the  opposite  siga  For  example, 

y,  =  -a  +  bi  y^  =  -a-bi 

y>*N=-b+ai  y,,^=+b+ai 

For  degenerate  modes,  the  real  parts  are  alwt^  zero,  and  the  coigugate  imaginaiy  parts  are 

indq>aident. 

y^  =  0+ci 

All  in  all,  there  are  only  two  unique  values  to  plot  for  each  pair  of  conjugate  modal 
variables.  For  die  non-degenerate  modes  we  plot  the  real  vs.  imaginaiy  parts  of  the  modal 
variables  to  get  phase  portraits.  Thoe  only  needs  to  be  one  plot  for  each  coiyugate  pair  due  to 
the  symmetiy  (i.e.  pick  y^  or  y^^^ . ).  Fix’  die  d^enerate  modes,  we  plot  the  imaginary  parts 
of  the  of  the  conjugate  y^ ,  y,^,^  variables.  In  the  linear  regime  of  the  periodic  orbit,  these 
phase  portraits  have  the  characteristic  forms  of  centers,  saddle  points,  and  straight  lines.  The 
centers  arise  fixim  modes  corresponding  to  purdy  imaginaiy  Poincar^  eiqxments,  the  saddles 
fixim  modes  with  real  Poincare  exponents,  and  the  straight  lines  fixxn  degenerate  modes  (zero 
Poincare  exponoits).  For  our  periodic  orbit  case,  there  are  two  modes  with  purely  imaginaiy 
Poincare  exponents  and  one  degenerate  motfo.  As  long  as  we  st^  in  the  linear  regime,  the 
three  plots  should  nominally  show  two  circles  and  one  line. 

The  behavior  of  die  three  modes  maps  back  to  classical  orbital  characteristics  in  the 
pt^ical  variables.  For  the  degenerate  mode  (mode  1),  variable  y^  correqxxids  to  a  time 
displacement  of  the  periodic  orbit,  while  y^  corresponds  to  negative  energy  (Wiesel, 

1981  ;236;  Wiesel,  1993:23).  We  refer  to  mode  1  as  the  time/emrgy  mode.  Modes  2  and  3  are 
die  vertical  and  planar  modes  respectively,  which  can  be  observed  1^  looking  at  the  eigenvector 
matrix  F  in  Appendix  A.  Vectors  2  and  S  have  only  vertical  components  (  z  and  p, ),  while 
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vectors  3  and  6  have  no  verticat  components.  The  Poincare  exponents  (or  modal  fi^uencies) 
of  the  these  two  modes,  "when  refo^enced  to  innlial  space,  imply  modal  periods  of  8.72S  yr  for 
the  planar  mode,  and  18.703  yr  for  the  vertical  mode.  These  are  the  fiuniliar  periods  for  the 
advance  of  the  perigee  and  regression  of  the  node  in  the  lunar  theory. "  (Wiesel,  1 98 1  :S88-589) 
Therefore,  we  call  mode  2  the  inclination  mode  and  mode  3  the  eccentricity  mode. 

To  demonstrate  the  characteristics  of  these  modal  plots,  we  start  1^  picking  tlw  near- 
periodic  orbit  trajectory  to  be  the  actual  periodic  orbit  we  found  earlier.  This  should  make 
Sx(t)  zero  for  all  time,  which  makes  the  modal  vector  zero  as  well.  The  modal  plots  are 
simply  points  with  random  noise  at  (0,0)  (see  Figures  4-6).  Next,  we  excite  the  modes  by 
perturbing  the  near-periodic  orbit,  x  ,  to  something  other  than  the  periodic  orbit.  This  is 
simply  a  matter  of  changing  die  initial  conditions.  For  illustration  purposes,  each  mode  is 
excited  independently  with  appropriate  initial  conditions.  Figures  7-9  show  examples  of  modal 
plots  from  orbit  tnyectories  in  the  linear  regime  of  the  periodic  orbit.  They  are  nearly  perfect 
circles  (Figures  8  and  9)  and  a  somediing  approaching  a  straight  line  (Figure  7)  when  the  initial 
conditions  are  perturbed  within  the  linear  region.  Figures  10-12  show  the  modes  of  an  orbit 
well  outside  die  linear  regime.  They  are  unstable  spirals  and  erratic  oscillations. 

In  the  next  section,  we  want  to  determine  if  die  actual  lunar  orbit  is  legitimately  one  of 
the  near-periodic  orbits.  In  other  words,  is  the  moon's  orbit  in  or  near  the  linear  r^me  of  our 
periodic  orbit? 
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Figure  4.  Time/Energy  Mode  using  Periodic  Orbit  Initial  Conditions. 
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Mode  3:  Not  Excited 


Figure  6.  Eccentricity  Mode  using  Periodic  Orbit  Initial  Conditions. 
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Figure  7.  Time/Energy  Mode  using  Initial  Conditions  in  the  Linear  Region. 
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Mode  2:  Linear  Regime 


Figure  8.  Inclination  Mode  using  Initial  Conditions  in  the  Linear  Region. 
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Figure  9.  Eccentricity  Mode  using  Initial  Conditions  in  the  Linear  Region. 


Mode  1:  Beyond  Linear  Region 


Figure  10.  Time/Energy  Mode  using  Initial  Conditions  Beyond  the  Linear  Region. 
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Figure  12.  Eccentricity  Mode  using  Initial  Conditions  Beyond  the  Linear  Region. 
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3.7  Actual  Lunar  Orbit 

We  now  want  the  near-periodic  orbit  to  be  the  actual  moon's  orbit  This  section 
discusses  the  method  for  using  real  historical  lunar  data  to  get  an  initial  condition  state  vector. 
jc(0) ,  for  modal  development  With  rrference  to  Pohlen  (1992),  this  development  is  what  he 
would  call  the  "exact  representation”  of  the  modal  variables. 


3.7.1  Getting  x(0)  Jromthe EjAtemeris 

To  begin,  we  select  some  random  time  at  which  a  new  moon  occurred,  and  take 
ephemerides  from  five  evenly  spaced  times  that  surround  the  new  moon  Specifically,  die 
following  data  is  collected  firom  the  American  EfAiemeris  and  Nautical  Almanac  (U.S.  Naval 
Observatoiy,  1%7)  at  each  of  the  five  times: 


lat^  •  Apparent  Latitude  of  die  Moon  with  respect  to  the  Ecliptic  Plane 

/on.  -  Apparent  Loi^tude  of  the  Moon  from  the  Vernal  Equinox  in  the  Ecliptic  Plane 

77r  ~  Horizontal  Parallax  of  the  Moon 

Ion,  -  Longitude  of  the  Sun  from  die  Vernal  Equinox  in  the  Ecliptic  Plane 


The  horizontal  parallax,  shown  in  Figure  13,  tmslates  to  a  distance,  R ,  between  the  earth  and 
moon  as 


R  = 


sinfl, 


(70) 


The  data  gathered  gives  information  that  looks  like  Figure  14,  but  what  is  really  needed 
is  the  position  of  the  moon  widi  respect  to  the  earth  in  a  coordinate  system  that  rotates  with  the 

A  A  A 

earth  around  the  sun.  Figure  IS  shows  sudi  a  geocentric  coordinate  system,  b  ,  vihose  b,  -bj 

A 

plane  sU^  in  the  ecliptic  and  whose  b,-axis  alw^  points  towards  the  sun.  This  coordinate 
frame  is  referred  to  as  the  Earth  Centered  Rotating  (ECR)  firame. 
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Earth 


T 

Fifac  Pbint  of  Ades 

Figure  14.  Data  Gadwred  fiom  Ephemerides. 
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Figure  1 S.  Earth  Centered  Rotating  (ECR)  Coordinate  Frame. 

Subtracting  the  sun's  longitude  fiom  the  moon's  longitude 

Slon^=lon^-lon,  (71) 

yields  the  desir'd  moon  position  information  with  respect  to  new  ECR  frame  in  the  polar 
coordinates:  lat^ , Slon^.R  (see  Figure  16).  To  g^  from  polar  coordinates  to  rectangular 
coordinates  we  use 

jc'  Rcos{latJ)co^SlonJ) 

=  .F'  =  Rco^kitJ^sin{SlonJ)  (72) 

z'  Rsin{latjj 

Figure  17  shows  this  transformation. 


39 


Figure  17.  Rectangular  ECR  Coordinates  fiom  the  Polar  Coordinates. 
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The  last  coordinate  transfixmatioa  is  a  translation  from  the  Earth  Centered  Rotating 
coordinate  frame  to  the  restricted  three  body  (R3B)  coordinate  frame.  This  is  accomplished 
simply  by  subtracting  die  distance  between  the  two  frames  from  the  ECR  position  vector 


^RiB 


0 

0 


(73) 


We  now  have  the  moon's  position,  fm  each  of  the  five  data  points  (times)  near  our 
chosen  new  moon,  in  die  restricted  three  body  problem  coordinate  systan.  A  good 
approximation  for  the  moon's  position  at  the  time  of  new  nxxxi,  r  (0) ,  is  derived  from  these 
five  surrounding  positions  using  the  Lagrange  Interpolation  Formulas  presented  in  the 
Handbook  of  Mathematical  Functions  (Abramowitz  and  Stegun,  1965:879).  For  our  five 
point  interpolation. 


'(o)  =  r(uf,  +ah) » 


(a*-l)a(a-2)  (a  -  l)o(a*  -  4)  _ 
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("-i) 


(a^-l)(a'-4)  (a  +  l)a(a'-4) 

+ - : - r{Uo) - 7 - r{u,) 


(a^-l)a(a  +  2) 

+ - - rM 


(74) 


vtdiere  h  is  die  stq)  size  between  each  of  the  five  ephoneris  times,  u,  are  the  five  equally 
spaced  ephemeiis  times,  r(u)  are  the  position  vectors  found  firom  £q  (73),  and  a  is  the  time 
difference  betweoi  the  third  data  point  ( Ug )  and  the  time  of  new  moon. 

To  get  die  initial  velocity  vector  at  the  new  moon  time  r  (0) ,  we  take  die  same  five 
data  points  and  use  die  derivative  form  of  Lagrange's  Formula  (Mramowitz  and  Stegun, 
1965:883) 
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F(o)  =  r(uo+ah)»- 


H 


2a* -3a* -a  +  l_/  v  4a* -3a* -8a +  4  _/  ^ 

- rly  i) - .) 

12  '  *'  6  ' 


2a*-5a_/  \  4a* +3a* -8a-4  _/  \ 
+ - 1 - 'M - 7 - '*(".) 

2  O 


+ 


2a*+3a*-a-l 


12 


(75) 


Since  r(/)  and  g(t)  are  the  same  thing,  the  initial  inertial  velocity  vector,  or  generalized 
momenta  vector,  p(0)  is  calculated  using  Eqs  (13)-( IS)  as 


p(0)  = 


x-y 
y  +  x 
z 


(76) 


Finally,  we  have  an  initial  state  vector  for  the  actual  moon 


(qm 


(77) 


U(o); 

This  state  vector  can  now  be  used  to  devek^  the  modal  variables  (§3.6)  and  determine  if  foe 
real  moon  is  in/near  the  linear  regime  of  our  p«iodic  orbit 


3.7.2  Preliminary  Results 

* 

It  should  not  matter  which  new  moon  we  select  from  foe  ephemerides,  therefim  we 
pick  one  at  random.  A  real  new  moon  occurred  on  9  Feb  1%7  at  1044  UT,  so  foe  five  data 
times  are  spaced  one  day  {4>art  starting  on  7  Fd>  at  0000  ITT.  Hie  interpolatfon  variables  are 
shown  below  and  foe  ^hemerides  are  summarized  in  ^jpendix  B. 

644.0 

fr=1.0  (days)  a  = -  (dttys) 

^  ^  '  1440.0 
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tt_j  =  7  Feb  OOOOUT 
=  8  Feb  OOOOUT 
Mo  =  9  Feb  OOOOUT 
M,  =10  Feb  OOOOUT 
Mj  =11  Feb  OOOOUT 


The  interpolatioa  method  returns  a  restricted  three  body  state  vector  for  the  new  nnoon 
that  looks  like 


-.997341 

-4.78964E-7 

-2.29153E-4 


jc(0)  = 


1.10990E-3 


-.966821 


3.79373E-4 


(77) 


After  doing  the  modal  development,  the  modal  variables  are  plotted  md  shown  in  Figures  18, 
19,  and  20.  At  best,  we  expect  to  see  circles  fiar  the  inclination  and  eccentricity  modes  and  a 
straight  line  for  the  time/energy  mode.  At  worst,  vn  expect  some  higher  firequency 
perturbations  about  the  two  circles  and  oscillating  (but  bounded)  time  and  energy. 
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Variable  y4 


Mode  1 :  Preliminary  Lunar  Orbit 


Figure  18.  Time/Energy  Mode  using  Preliminary  Lunar  Orbit  Initial  Conditions. 
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Figure  20.  Eccentricity  Mode  using  Preliminary  Lunar  Orbit  Initial  Conditions. 
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Obviously,  there  is  a  problem.  While  mode  3  looks  somev^tat  close  to  a  circle  with 
perturbations,  mode  2  spirals  completely  out  of  the  periodic  orbit’s  linear  r^on.  In  addition, 
the  energy  (variable  >'4)  of  the  system  seems  to  be  going  unbounded.  To  check  these  results, 
we  choose  another  new  moon  from  the  ephemerides;  9  Feb  1986.  These  efhemerides, 
interpolation  variables,  and  initial  state  vector  are  in  Appendix  C.  The  results  from  die  Feb  86 
data  are  veiy  similar  to  the  Feb  67  results. 

The  problem  is  two  fold;  1)  we  are  using  initial  conditions  derived  from  real  data  for  a 
particular  month,  but  the  pericxl  we  are  using  is  the  average  lunar  period  over  histoiy.  The 
actual  lunar  period  during  Feb  67  is  longer  than  the  average  value;  2)  for  die  initial  state  vector 


weknowdiat  ^(0)  correctly  matches  the  real  world,  but  since  p(0)  is  derived  ^-om  the  real 
world  and  our  R3B  model  doesn't  include  eveiything  die  real  world  does,  then  p  (0)  does  not 
provide  our  model  widi  the  initial  conditions  necessary  to  r^roduce  the  moon's  orbit. 

To  better  illustrate  the  errors,  reference  Figure  21 .  This  shows  the  X  and  Y 
components  of  the  pericxlic  orbit  and  die  real  orbit  in  physical,  R3B  ccxirdinates.  Both  were 
integrated  using  the  R3B  model  for  the  same  amount  of  time.  Th^  are  siqiposed  to  have  the 
same  period,  but  i^e  die  pericxlic  orbit  makes  it  back  to  the  X-axis,  die  real  orbit  does  not. 
Our  dynamic  mcxlel  and  average  pericxl  are  conflicting  with  what  the  ^hemeris  says  the  real 
orbit  should  be.  Therefore  the  problem  is  this;  givm  that  q  {Q  is  correct,  find  piQ  such 
that  our  R3B  model  will  rqxoduce  subsequent  q(tys  that  match  real  data. 
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Figure  21.  Periodic  Orbit  and  Preliminary  Lunar  Orbit  in  Physical  Space.  X  and 
Y  components  of  the  Restricted  Three  Body  Problem. 


3. 7. 3  Fitting  the  Initial  Momenta  to  Our  Dynamics  Model 

The  solution  to  the  problem  above  is  to  iteratively  fit  the  initial  momenta  such  that  die 
integrated  position  matches  the  desired  real  position.  The  process  is  similar  to  finding  the 
periodic  orbit  initial  conditions  (§3.2). 

We  start  by  picking  some  later  time,  r, ,  after  the  new  moon  and  cdlecting  more  data. 
The  ^hemeris  is  converted  with  the  same  coordinate  transfwmations  explained  in  section 
3.7. 1  to  get  the  position  vector  ^(/,) .  The  initial  new  moon  state  vector  (Eq  (77))  is 
integrated  forward  in  time  to  /,  so  that  the  error  in  position  at  /,  is  calculated. 


Now  Eq  (27)  is  again  recalled. 


Sx{t,)  =  <l>{t^,tJSx{to) 


(79) 


MiiV 

^pOi)^ 

O. 


nj 


(80) 


so  if  we  assume  that  the  initial  position  is  correct  (  5r{tg )  =  0  )  and  the  dififerences  in  final 
momenta  (  Sp{t^ ) )  are  uninqx)rtant,  then  we  can  isolate  only  the  pails  that  concern  us  to  get 


=  sp{0 

-«(/,)  =  0,2  Sp{t^) 

Solving  for  Sp(tg)  yields  die  correction  to  die  initial  momenta  vector  so  that 

pOo)^  =  pUo)m-^^p(^o) 


(81) 

(82) 


(83) 


This  new  initial  momenta  vector  is  substituted  back  into  the  initial  state  vector  and  the  cycle  is 
rqjeated  until  die  error  meets  some  tolerance  (we  used  10-'°).  When  the  error  is  small 
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enough,  we  have  a  set  of  initial  conditions  that,  when  int^rated  forward  with  our  R3B 
dynamics  model,  echo  the  actual  lunar  qrfiemeris  for  time  t, .  We  then  pick  a  new  value  for  t, 
,  and  use  those  initial  conditions  as  a  new  first  guess.  The  process  is  continued  imtil  a  large 
enough  t,  is  found  such  that  the  average  period  of  the  real  orbit  matches  the  period  of  the 
periodic  orbit. 

The  end  product  here  is  an  initial  state  vector  that  will  yidd  the  exact  modal 
representation  of  the  lunar  orbh. 
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3.8  Perturbation  Theory 


In  the  last  sectkMi,  we  developed  the  exact  modal  representation  ci  the  moon's  orbit. 
Here,  the  goal  is  to  develop  an  expanded  modal  solutioa  That  is,  we  want  the  equations  of 
motion  for  the  modal  variables  based  on  the  perbirbation  expansion  of  the  new  modal 
Hamiltonian. 

We  start  with  the  original  Hamiltonian,  X(jc)  ,  and  expand  it  in  a  Tqrlor’s  swes  about 
the  periodic  orbit  (i.e.  x  =  x^  or  Sx  =  Q  ). 

»(*)=*. +«,+»,+  - 


d% 


dx.. 


Sx.. 


Si=0 


I  «  6 


2!  dx^dx^ 


Sx^Sx^ 


Sx^ 


1  «  6  « 
«.=izzz 


3!  dx^dx^dx^ 


Sx  Sx^Sx^ 


6x^ 


(84) 


The  tensor  notation  for  Eq  (84)  is 


JC{x)  =  Xo  +X^^  Sx^  Sx^^X^SXj^  +  •••  (85) 

The  first  term,  ,  is  a  scalar  constant  because  it  is  the  Hamiltonian  fiv  the  pmodic  orbit 

For  simplicity,  we  choose  the  scalar  constant  to  be  zero.  'The  second,  or  linear  term  is 
identically  zero,  because  it  describes  the  motion  of  the  periodic  trajectory  with  respect  to  itself 
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The  third,  or  quadratic  term  is  the  Floquet  problem,  and  becomes  a  constant  co^cient,  linear 
systrai  in  the  new  variables"  (Ross,  1991:31).  Therrfore,  the  expanded  Hamiltonian  now 
looks  like 

+  ^X^^Sx^Sx^Sxj^Sx^  +  •  (86) 

By  writing  the  transformation  to  modal  variables  Sx  =  Ey  in  tensor  notation  Sx^  -^^yj  > 
we  can  canonically  transfixm  Eq  (86)  to  die  new  expanded  Hamiltonian  in  modal  variables. 


X{y)=j^X^y,yj  +^X^y,y,y,  +^*4*.  y,y,y,y,  +-  (»T) 

i^iere 

Xu==X^,e^e,  (8S) 

(89) 

The  term  is  the  constant,  symmetric  matrix  S  introduced  earlier  (§3.3).  The  ^3^  , 

^4^  ,  and  higher  terms  are  periodic  tensors  since  the  partials  tensors  ( X  terms)  and  the  E 
matrix  are  both  periodic  with  the  same  period  as  our  periodic  orbit  We  get  the  pardals  tensors 
as  well  as  the  elements  of  the  E  matrix  fiom  the  Fourier  coefficients  of  x^{t)  and  E{t)  just  as 
befixe. 

As  mentioned  in  section  3.4,  making  y{t)  =  E~'{t)Sx{t)  a  canonical  transformation 
ensures  that  die  new  modal  Hamiltonian  X{y)  fiollows  Hamilton's  equations.  This  is  useful, 
because  we  want  modal  equations  of  motion  and  we  get  them  fiom  Hanuhon's  equations; 
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Again,  using  tensor  notation,  we  get 

The  first  term  is  the  constant  coefficient,  linear  system  we  estid>lished  with  Eq  (48),  where 
=  ZS  =  J.  By  itself^  ~  should  create  modal  plots  that  are  perfect  circles, 

saddles,  and  straight  lines  fix  the  respective  modes.  The  second,  third,  and  subsequent  terms 
represent  the  perturbations  off  of  the  linear  systera 

We  now  have  equations  of  motion  in  tmas  of  the  modal  variables,  with  initial 
conditions  derived  1^  canonically  transfixming  the  initial  state  vector  fimn  physical  variables. 
The  modal  plots  fiom  this  expansion  solution  should  compare  closely  with  those  firom  the  exact 
solution  discussed  in  the  last  section. 

Specifically  in  our  case,  when  the  e}q)ansion  only  includes  the  Xj  term,  we  should  get 
perfect  circles  for  the  inclirudion  and  eccentricity  modes  aiKi  a  straight  line  fix  the  time/energy 
mode.  As  each  of  die  additional  poturbation  terms  are  included  in  die  expansion,  the  modal 
plots  should  agree  more  and  more  dosdy.  Ifallofdie  most  significant  terms  are  included,  the 
two  modal  representations  should  match  well.  As  Fouria-  series  computations  get  extremely 
cumbersome  beyond  the  third  term,  we  choose  to  limit  our  expansion  to  . 


IV.  Software 


To  do  this  study,  several  FORTRAN  and  MATLAB™  computer  codes  were  used. 

This  chapter  discusses  the  specific  software  that  was  (teveloped  to)  accomplish  each  of  the 
corresponding  sections  of  Chapter  Three.  Aldiough  the  software  codes  are  not  presorted  in 
this  document,  they  can  be  r^ioduced  from  the  information  in  Chapter  Three  or  obtained  from 
die  author  or  Dr.  Wiesel.  Unless  otherwise  stated,  each  code  was  written  in  FORTRAN  77. 

4.1  Periodic  Orbit 

The  program  PERIOD  takes  our  initial  conditioos  guess  frmn  section  3.2  (Eq  (37)) 
and,  using  a  predictor/corrector  algorithm  called  HAMING,  iterates  toward  a  nearby  set  of 
initial  conditions  that  whoi  integrated  for  one  period,  creates  a  periodic  orbit  PERIOD  is  set 
up  specifically  to  handle  symmetric  periodic  orbits  in  the  restricted  three  botfy  problem.  The 
iteration  process  was  discussed  in  more  detail  in  section  3.2.2.  After  finding  dw  periodic  orbit 
initial  conditions,  PERIOD  calculates  die  moiKxlron^  matrix,  eigenvalues,  eigenvectors  ( F 
matrix),  and  Poincare  exponents  ( J  matrix)  as  described  in  section  3.3. 

4.2  Symplectic  Normalization 

Before  doing  the  symplectic  normalization,  the  RPEXVEC  code  is  needed  to  determine 
die  restricted  problem  extended  (or  generalized)  eigenvector  for  those  systems,  like  ours,  with 
degoierate  modes.  RPEXVEC  uses  hardwired  values  firom  selected  rows  of  die  monodimny 
matrix  ( O )  and  die  conjugate  eigenvector  (  /,  )  to  solve  Eq  (67)  for  the  extended  eigenvector 
/;  using  a  linear  equations  solver.  The  eigenvector  matrix  F  is  modified  accordingly  and 
fed  into  program  SYMNRM. 

SYMNRM  does  the  actual  symplectic  normalization  of  F  to  E  by  calculating  the 
diagonal  entries  of  the  D  matrix  (Eq  (63)).  Hie  ouqiut  of  SYMNRM  is  the  initial  E  matrix. 
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4.3  Modal  Variable  Development 

The  periodic  orbit  initial  conditions,  the  J  matrix,  and  the  initial  E  matrix  are  all 
inputs  to  the  MOOEDV  pn^ram.  MODEDV  integrates  the  periodic  orbit  and  £  matrix 
forward  in  time  to  one  period,  saving  values  of  each  at  evenly  spaced  intervals.  Using  the 
harmonic  analysis  subroutine  CMPFOR  (based  on  the  algorithm  of  Brouwer  and  Clemence 
(1961;108-1 13)),  we  then  get  a  set  of  coii4)lex  ^qxxiential  Fourier  coefficients  for  both  the 
periodic  orbit  state  vector  and  £  (§  3.6). 

The  Fourier  coefficients  are  loaded  into  the  program  TRAJ,  as  is  the  initial  condition 
state  vector  of  some  near-periodic  tnyectory.  TRAJ  integrates  the  near-periodic  tnyectory  in 
discrete  steps.  At  each  step,  TRAJ  calls  the  subroutine  MODEVR  which  sums  the  Fourier 
coefficients  of  the  periodic  orbit  and  E  ,  and  calculates  the  modal  vector  using  Eq  (69).  The 
result  is  a  series  of  modal  vectors  evenly  spaced  throughout  the  integr^on  time.  These  modal 
vectors  are  dien  easily  plotted  as  discussed  in  section  3.6. 

4. 4  Real  Moon  Initial  Conditions  from  Data  (Exact  Solution) 

A  M ATLAB™  code  called  FEB67.m  is  used  to  calculate  the  real  moon's  initial 
conditions  from  February  1967  data  (an  identical  code  called  FEB86.m  was  also  created  for 
1986  data).  The  sun/moon  information  is  loaded  for  five  evenly  spaced  times  around  the  new 
moon,  transformed  into  the  Earth  Coitered  Rotating  frame,  and  finally  translated  to  the 
restricted  dvee  body  frame.  FEB67.m  dioi  interpolates  to  get  the  moon's  position  and  velocity 
vector's  as  shown  in  section  3.7.1.  Afrer  doing  one  last  transformatirxi  of  die  velocity  vector  to 
the  inertial  velocity  vector,  FEB67.m  returns  an  initial  condition  state  vector  for  the  actual 
moon.  This  initial  condition  state  vector  is  nominally  loaded  into  program  TRAJ  as  the  near- 
periodic  trajectory.  However,  it  was  shown  in  section  3.7.2  that  the  initial  momenta  needs  to 
be  fitted  to  our  dynamics  model. 
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Program  FTTICS  takes  dw  initial  new  moon  state  vector  (output  from  FEB67.m),  and 
iterates  toward  an  initial  momenta  vector  such  that  the  integrated  position  matches  the  position 
listed  in  the  ephemerides  at  some  time  after  new  nvxxi.  The  process  is  similar  to  that  used  for 
program  PERIOD  and  the  details  are  shown  in  section  3.7.3.  The  output  of  FITICS  is  the  new 
initial  condition  state  vector  that  is  entered  into  TRAJ  as  the  near-periodic  trajectcwy. 

4.5  Expanded  Modal  Solution 

There  are  three  main  programs  used  to  create  the  modal  expansion  plots;  K3DEV, 
K4DEV,  and  YTRAJ.  K3DEV  develops  the  third  order  periodic  tensor  ( )  of  the 
expanded  Hamiltonian  in  modal  variables.  At  eadr  evenly  spaced  time,  K3DEV  reads  in  and 
sums  the  Fourier  coefficients  ft>r  the  periodic  orbit  state  and  E .  Then,  the  third  order  partial 
tensor  of  the  original  Hamiltonian  ( )  is  retrieved  from  subroutine  ORDER3.  Using  the 
summation  convention  shown  in  section  3.8,  we  get  the  third  order  tensor  ( ^3^  )  at  evenly 
spaced  times  duoughout  the  period.  Again  using  the  harmonic  analysis  subroutine  CMPFOR, 
we  get  a  series  of  Fourier  coefficioits  for  ^3^  . 

K4DEV  does  the  same  thing  as  K3DEV  except  that  it  uses  subroutine  ORDER4  and 
finds  the  faurih  order  periodic  tensor,  ^4^  ,  of  the  expanded  Hamiltonian  in  nradal  variables. 

Program  YTRAJ  is  the  counterpart  of  TRAJ.  Instead  of  integrating  the  orbit  in 
physical  variables  and  transforming  to  modal  variables  at  each  stq),  YTRAJ  integrates  the 
rtKxlal  variables  directly.  The  initial  state  vector  is  transformed  from  physical  variables  to 
modal  variables,  and  the  ^3  /  Fourier  coefficients  are  read  in.  At  each  time  step,  the 
and  X4  Fourier  coefficients  are  summed  and  the  modal  variables  are  integrated  by  Eq  (92). 
The  result  of  YTRAJ  is  a  series  of  modal  vectors  that,  when  plotted  over  time,  should  match 
the  plotted  res  'Its  of  TRAJ. 
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V.  Results  and  Discussion 


5.1  Exact  Modal  Results 

In  section  3.7.1,  we  find  an  initial  state  vector  in  physical  variables  fw  the  real  moon, 
however,  the  initial  momoita  needs  to  be  fitted  such  that  the  modd  integrations  match  the  real 
data  at  a  future  time  (section  3.7.3).  To  do  die  fitting,  values  of  t,  are  selected  starting  widi 
two  weeks  after  the  new  nKXMt  Next,  one  month  is  tried,  then  six  months,  and  finally  one  year 
after  die  new  moon.  The  result  is  a  well  approximated  initial  condition  state  vector  in  {diysical 
variables  for  the  moon.  This  inidai  lunar  state  vector  is  the  near-periodic  orbit  that  is  integrated 
and  transformed  to  modal  variables  using  the  TRAJ  program  The  modes  produced  are  what 
we  call  the  exact  nwdal  representation. 

As  stated  earlier,  die  first  result  of  interest  is  whether  the  lunar  orbit  is  in  or  near  the 
linear  region  of  the  periodic  orbit  The  modal  pk>ts  of  the  exact  modal  representation  are 
shown  in  Figures  22, 23,  and  24.  Figure  22  is  a  plot  of  the  tiine/energy  mode  (nnode  1).  The 
epoch  time  displacement  (variable  .y, )  is  shown  oscillating  back  and  fiarth,  which  basically 
represents  in-track  fluctuations  of  die  morxi's  orbit  about  the  original  pei.odic  orbit  The 
system  energy  (variable  .^4 )  is  also  shown  oscillating,  however  it  is  bounded.  Of  the  three 
modes,  mode  1  has  the  least  stringent  requiremoits  for  acceptable  bdiavior.  As  a  general  rule, 
according  to  Dr.  Wiesel,  the  time/energy  mode  is  considered  to  be  behaving  reasonably  as  long 
as  the  energy  stays  bounded.  As  you  can  see,  it  does  stays  bounded. 

The  next  plot  is  the  inclination  mode,  or  mode  2,  shown  in  Figure  23.  We  are  looking 
for  something  close  to  a  circle,  and  that  is  basically  what  we  get  The  plot  shows  higher 
ftequmcy  oscilladrms  about  a  generally  circular  modoa  These  hi^ier  frequency  loops 
characterize  die  perturbations  to  the  linear  Floquet  periodic  orbit.  In  fact,  the  number  of  higher 
frequency  loops  in  the  pattern  (three)  correspoids  to  the  order  of  dominant  terms  in  the  Fourier 
series.  Chven  the  relative  size  of  the  oscillations,  we  can  conclude  diat  this  mode  is  not 
necessarily  in  the  linear  regime,  but  it  is  definitely  close. 
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Figure  22.  Exact  Modal  Representation  of  the  Time/Energy  Mode  using  Fitted 
Initial  Conditions. 
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Mode  2:  Exact  Modal  Representation 


Figure  23.  Exact  Modal  Representation  of  the  Inclination  Mode  using  Fitted  Initial 
Conditions. 
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Figure  24.  Exact  Modal  Representation  of  the  Eccentricity  Mode  using  Fitted  Ini¬ 
tial  Conditions. 
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Figure  24  shows  the  plot  for  mode  3,  the  eccentricity  mode.  As  with  mode  2,  we 
e}q)ect  to  see  a  circle  and  again,  the  plot  shows  a  circle  widi  higher  frequency  oscillations.  In 
this  case,  however,  the  eccoitricity  mode's  motion  is  definitely  closer  to  a  perfect  circle. 
Therefore,  mode  3  is  close  to  the  linear  regime,  but  not  quite  in  it. 

Looking  at  all  three  modal  plots,  it  is  concluded  that  the  moon's  real  orbit  is  not  in  the 
linear  region  of  foe  periodic  orbit,  but  it  is  close  «K>ugh  that  we  can  see  the  reference  linear 
behavior  about  which  foe  higher  frequency  oscillations  are  being  perturbed.  The  next  step  is  to 
proceed  with  foe  rest  of  foe  analysis  and  compare  these  exact  nmxlal  plots  with  foe  expanded 
modal  plots.  They  should  match  closely. 

5.2  Expanded  Modal  Results 

In  this  section  we  present  foe  modal  plots  that  are  derived  firom  foe  modal  equations  of 
motion  based  on  foe  perturbation  expansion  of  foe  new  modal  Hamiltonian  (Eq  (92)).  These 
plots  are  foe  results  of  program  YTRAJ,  and  we  call  them  foe  expanded  modal  rq}resentation. 
They  are  compared  with  foe  exact  modal  representatitm  results  from  foe  last  section.  The  two 
nKxlal  representations  should  match  closely  if  we  include  enough  perturbation  expansion  terms. 
The  comparison  accomplishes  three  things. 

1 .  It  provides  redundancy  of  the  exact  modal  plots  to  satisfy  ourselves 
that  foe  technique  was  properly  applied. 

2.  It  verifies  foat  we  have  successfully  transformed  the  Hamiltonian 
itself  to  modal  coordinates. 

3 .  It  determines  if  foe  perturbatitms  are  mostly  (or  entirely)  derived 
from  only  foe  X  -terms  that  we  choose  to  include  in  foe  expansion. 

We  first  run  foe  YTRAJ  program  only  considering  foe  first  term  ( )  of  foe 
expanded  representation.  The  expanded  modal  plots  are  overlaid  with  foe  exact  modal  plots 


61 


and  die  results  are  shown  in  Figures  2S,  26,  and  27.  As  discussed  in  section  3.8,  the  X2 
expansion  term  by  itself  represents  a  constant  coefficioit,  linear  system.  It  is  the  Floquet 
portion  of  the  expansion  and  it  exhibits  perfect  linear  behavior.  For  modes  2  and  3  (Figures  26 
and  27),  we  get  a  perfect  circle  that  oveiit^  nicely  with  the  exact  modal  plot  to  show  die  path 
about  which  diey  are  perturbed.  Mode  1  (Figure  27)  is  a  perfectly  straight  line  as  expected, 
where  energy  is  fixed.  If  the  lunar  orbit  was  in  die  linear  regime  of  the  periodic  orbit,  the  X2 
term  would  essentially  be  the  only  term  in  the  expansion  necessaiy  to  match  the  two 
representations.  But  since  die  moon  is  definitely  outside  the  linear  region,  we  need  to  include 
the  perturbation  terms.  As  we  include  more  term  in  the  expansion,  we  should  start  perturbing 
off  of  the  perfect  circles/line,  and  start  matching  the  exact  modal  plots  more  closely. 

The  YTRAJ  program  is  run  again,  but  this  time  the  X^  term  is  added  to  the 
expansion.  These  modal  plots  are  overlaid  wndi  the  exact  modal  plots  and  shown  in  Figures 
28, 29,  and  30.  The  results  are  unreascxiable  in  that  the  expanded  plots  initially  follow  the 
exact  solution  plots,  but  then  spiral  off  toward  infinity.  For  example.  Figure  29  shows  diat 
mode  2  stiQrs  widi  the  first  high  fi-equency  perturbation  loop  for  a  while,  but  then  quickly 
spirals  away.  The  best  thing  to  do  is  to  ccmtinue  adding  perturbation  terms  to  the  expansion. 

Now  YTRAJ  is  run  widi  all  terms  in  the  e}q)ansion  out  to  the  X^  term.  The  modal 
plots  are  again  overlaid  with  the  exact  solution  and  presented  in  Figures  3 1, 32,  and  33.  They 
show  definite  improvement  over  the  previous  two  sets  of  modal  plots.  The  e}q>ansion  modes 
stay  with  the  exact  modes  a  lot  longer  this  time,  but  tiiey  still  eventually  spiral  off  toward 
infinity.  Again  using  mode  2  (Figure  32)  as  an  example,  die  behavior  matches  quite  well  for 
the  first  three  or  four  perturbation  loc^s  before  it  spirals  away.  The  conclusion  again,  is  that 
more  perturbation  terms  need  to  be  considered  in  die  expansion.  However,  instead  of  dealing 
with  a  computational  nightmare,  there  is  an  attractive  alternative. 

It  is  our  contention  diat  the  plots  would  match  if  we  added  enough  terms  to  die 
expansion,  but  because  the  real  moon  is  not  in  die  linear  region  of  the  periodic  orbit  it  takes  a 
great  number  of  the  expansion  terms  to  get  a  perfect  match.  Consider  diat  X^^  isa  6x6x6 


62 


Mode  1:Exact(Solid)  VS  Expanded  to  K2(Dotl6d) 
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Figure  25.  Exact  vs  Expanded  Modal  Solution  (K2  Only).  Time/Energy  Mode. 


Figure  26.  Exaet  vs  Expanded  Modal  Solution  (K2  Only).  Inclination  Mode. 
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Figure  27.  Exact  vs  Expanded  Modal  Solution  (K2  Only).  Eccentricity  Mode. 
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Figure  28.  Exact  vs  Expanded  Modal  Solution  (K2  &  K3  Only).  Time/Energy 
Mode. 
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Figure  29.  Exact  vs  Expanded  Modal  Solution  (K2  &  K3  Only).  Inclination  Mode. 
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Figure  30.  Exact  vs  Expanded  Modal  Solution  (K2  &  K3  Only).  Eccentricity 
Mode. 
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Figure  31.  Exact  vs  Expanded  Modal  Solution  (K2,  K3  k  K4).  Time/Energy 
Mode. 
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Figure  32.  Exact  vs  Expanded  Modal  Solution  (K2,  K3  &  K4).  Inclination  Mode. 
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Figure  33.  Exact  vs  Expanded  Modal  Solution  (K2,  K3  &  K4).  Eccentricity  Mode. 


tensor,  4  is  a  6x6x6x6  tensor,  and  that  each  element  of  those  t«isors  has  4 1  Fourier 

coefficients.  The  computational  requirements  of  the  6x6x6x6x6x41  (or  higher)  array 
needed  to  facilitate  the  ^3  (or  higher)  term  would  obviously  be  great. 

Rather  than  deal  with  the  laborious  task  of  including  die  X, ,  X^,  etc.  terms  to  get  the 
modal  plots  to  match,  we  created  a  fictitious  mocm  that  exhibits  the  same  basic  dynamics  of  die 
real  moon  but  is  much  closer  to  the  linear  regime  of  the  periodic  orbit.  To  get  this  fictitious 
moon,  we  decrease  the  amplitude  of  the  real  moon's  modes  by  taking  the  initial  conditions  in 
modal  variables  and  dividing  by  an  arbitrary  constant  (we  use  1 0.0 ).  The  initial  conditions  in 
physical  variables  are  calculated  from  Eq  (46)  as 

Jjc(0)  =  F(0)y(0)  (93) 

and  loaded  into  program  TRAJ.  The  result  is  a  set  of  exact  modal  plots  for  the  fictitious  moon. 
Meanwhile,  the  initial  conditions  in  modal  variables  are  loaded  into  program  YTRAJ,  where 
the  only  expansion  terms  included  are  X2  ,  X^,  and  ^4  .  The  result  is  a  set  of  expanded 
modal  plots  for  the  fictitious  moon.  These  new  modal  plots  are  compared  and  shown  in 
Figures  34,  35,  and  36. 

Obviously,  the  fictitious  moon  is  much  closer  than  the  real  moon  to  the  linear  region  of 
the  periodic  orbit.  In  fact,  mode  3  (Figure  36)  is  basically  a  perfect  circle.  Since  we  are  so 
much  closer  to  the  linear  regime,  we  expect  that  less  expansion  terms  are  required  to  get  the 
exact  and  expanded  results  to  match.  That  is  definitely  the  case  here.  The  exact  and  expanded 
modes  match  very  closely.  Mode  1  (Figure  34)  breaks  down  the  earliest  of  the  three  modes, 
but  even  so,  the  form  compares  close  enough  to  make  the  point  and  it  shows  behavior  that  is 
much  closer  to  pure  time  drift  with  bounded  oiergy  (straight  line).  Mode  2  matches  quite  well 
and  mode  3  is  essentially  a  perfect  comparison. 

One  brief  side  note  is  in  order.  Notice  that  even  with  the  fictitious  moon,  the 
perturbations  for  the  inclination  mode  and  the  system  energy  {y^  are  growing  with  time.  The 
most  likely  reason  for  this  is  that  a  value  larger  than  one  year  may  be  needed  for  /,  (§  S.  1 ). 
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Variable  y4 


Figure  34.  Exact  vs  Expanded  Modal  Solution  (K2,  K3  &  K4)  for  Fictitious  Moon. 
Time/Energy  Mode. 


Figure  35.  Exact  vs  Expanded  Modal  Solution  (K2,  K3  &  K4)  for  Fictitious  Moon. 
Inclination  Mode. 


Figure  36.  Exact  vs  Expanded  Modal  Solution  (K2,  K3  &  K4)  for  Fictitious  Moon. 
Eccentricity  Mode. 


Since  die  new  modal  plots  matched  so  well,  we  have  shown  that  the  canonical  Floqu^ 
perturbation  theory  technique  was  successfully  applied  to  the  moon's  orbit  and  that  Eq  (92)  is  a 
modal  solution  to  the  real  moon,  however  a  great  number  of  terms  would  be  required  in  die 
expansion. 

5.3  Complex  Conjugacy 

In  the  paper  by  Wiesel  and  Pohlen  (1992:9),  it  was  assumed  that  the  eigenvector 
matrix  of  a  Hamiltonian  system  must  always  contain  complex  conjugate  eigenvectors  to 
perform  a  canonical  transformation.  As  it  turns  out,  the  eigenvectors  of  the  moon's  monodromy 
matrix  does  not  contain  complex  conjugate  eigravectors.  It  is  concluded,  therefore,  that 
symplectic  normalization  does  not  necessarily  need  to  involve  maintaining  complex  conjugacy 
of  the  eigenvector's  from  F  to  E. 
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VI.  Conclusions  and  Recommendations 


The  main  conclusion  of  this  work  is  that  the  canonical  Floquet  dieory  technique  was 
successfully  applied  to  the  moon's  orbit.  This  study  showed  that  caixMiical  Floquet  perturbation 
theory  is  applicable  to  a  system  which,  besides  being  dynamically  more  interesting  than 
previous  applications,  is  also  constrained  to  obey  real  historical  data.  We  used  the  restricted 
three  body  problem  to  get  Hill's  lunar  orbit  as  a  stable  periodic  reference  orbit,  and  appUed 
canonical  Floquet  theory  using  interpolated  inidai  conditions  from  real  moon  ephemerides  of 
February  1967.  The  moon's  interpolated  initial  position  proved  to  be  accurate,  but  the 
interpolated  initial  momenta  had  to  be  modified  to  accommodate  the  differences  in  our 
restricted  dvee  body  model  and  the  real  world  dyiuimics  We  successfully  fit  the  initial 
momenta  such  diat  the  moon's  position,  when  integrated  forward  fi>r  one  year,  still  matched  the 

data 

The  modal  plots  of  the  real  moon's  orbit  showed  that  it  was  outside  the  linear  region  of 
Hill's  periodic  orbit,  but  still  close  enough  to  demonstrate  the  technique.  The  time/energy 
mode  oscillated,  but  the  energy  remained  bounded,  while  the  inclination  and  eccentricity  modes 
oscillated  with  highe'  frequency  perturbations  about  an  expected  perfect  circle.  We  derived  the 
modal  solution  using  both  an  exact  and  an  expanded  modal  representation.  The  Hamiltonian 
was  successfully  transformed  to  modal  variables  and  the  modal  equations  of  motion  were 
found  out  to  diree  expansion  terms.  We  truncated  the  expanded  solution  after  each  of  the  first 
three  terms  and  compared  the  modal  plots.  The  results  demonstrated  that  the  exact  and 
expanded  representations  will  agree  given  that  enough  expansion  terms  are  included.  For  the 
fictitious  moon,  three  terms  were  sufficient  to  consider  the  modal  EOM  a  viable  modal 
solution,  however  the  real  moon,  being  further  away  from  the  linear  region,  requires  many 
more  expansion  terms. 

In  future  studies,  the  modal  solution  found  here  can  be  used  as  the  reference  tr^ectoiy 
for  fiirdier  analytic  perturbation  applicaticHis  sudi  as  Von  Ziepel's  medKxl.  In  addition,  the 
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higher  order  expansion  terms  could  be  included  to  con^)l^ely  validate  our  contention.  Also, 
the  fitting  of  initial  mcHnenta  to  future  ephemerides  could  be  extended  b^rond  one  year.  It 
might  decrease  the  magnitude  of  the  higher  frequefK:y  perturbations. 
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Appendix  A:  Monodromy  Matrix,  Initial 
Eigenvector  Matrix,  and  J  Matrix 

The  foUowing  matrices,  used  for  Ftoqi^  theory,  help  describe  our  periodic  orbit  and 
were  output  from  the  prc^ram  PERIOD. 


The  moiKxiromy  matrix  <1>(  r.O); 
row  column 

1  1  1.095712688474569 

2  1  -20.370277922605400 

3  1  O.OOOOOOOOOOOOOOOE-HXX) 

4  1  274.031141304028700 

5  1  -1.269564310370324 

6  1  O.OOOOOOOOOOOOOOOE-KXX) 

1  2  6.294761 1 105 16258E-001 

2  2  8.048342665042469E-001 

3  2  O.OOOOOOOOOOOOOOOE-HXX) 

4  2  1.269564305657530 

5  2  -8.349576338746477 

6  2  O.OOOOOOOOOOOOOOOE-HXX) 

1  3  O.OOOOOOOOOOOOOOOE+OOO 

2  3  O.OOOOOOOOOOOOOOOE-HXX) 

3  3  8.601930291457586E-001 

4  3  O.OOOOOOOOOOOOOOOE-HXX) 

5  3  O.OOOOOOOOOOOOOOOE-HXX) 

6  3  -7.345815692997451 

1  4  4.745632093038844E-002 

2  4  -1.471358089142643E-002 

3  4  O.OOOOOOOOOOOOOOOE-HXX) 

4  4  1.095712688719114 

5  4  -6.294761 110578383E-001 

6  4  O.OOOOOOOOOOOOOOOE-HXX) 

1  5  1.471358087321641E-002 

2  5  -1.491295028851029 

3  5  O.OOOOOOOOOOOOOOOE+OOO 

4  5  20.370277922636880 

5  5  8.048342661624743E-001 

6  5  O.OOOOOOOOOOOOOOOE-HXX) 
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monodronv  matrix  (ctMit) 


row  column 

1  6  O.OOOOOOOOOOOOOOOE-KXX) 

2  6  O.OOOOOOOOOOOOOOOE-KXX) 

3  6  3.540354992250493E-002 

4  6  O.OOOOOOOOOOOOOOOE-KXX) 

5  6  O.OOOOOOOOOOOOOOOE-KXX) 

6  6  8.601930291462128E-001 


The  initial  eigenvector  matrix  F(o)  : 

row  column  ( real ,  imaginary ) 

1  1  (.8.192793633048179E-014^.007180267772390E-027) 

2  1  (-7.539018574763626E-002,1.667302866291538E-017) 

3  1  (-1.217710722943405E-019,4.972171594093748E-023) 

4  1  (l.OOOOOOOOOOOOOOO.O.OOOOOOOOOOOOOOOE-KXX)) 

5  1  (-1.361796347709330E-011,.1.305682296333588E-026) 

6  1  (-9.189976165277692E-019,2.610493483499758E-023) 

1  2  (-7.539018574755219E-002,-1.439135915971290E-012) 

2  2  (1.405S9021S26S010E-Oil.-1.021808214943804EK)01) 

3  2  (-3.808102176268219E-017, 2.61771 1356936501E-018) 

4  2  (-1.757523295343237E-010,6.646921145391658E-001) 

5  2  (l.OOOOOOOOOOOOOOO.O.OOOOOOOOOOOOOOOE-KXX)) 

6  2  (-1.99933283 1133622E-017,-4.141361260090141E-017) 

1  3  (8.609626586666082E-017,-1.470913273776291E-017) 

2  3  (-9.403731537416513E-017,-4.199563329593159E-016) 

3  3  (3.092693812092022E-014,-6.942300233683507E-002) 

4  3  (1.391287225318998E-015,6.660459595692267E-015) 

5  3  (-1.332598592779162E-015,3.703060854033629E-016) 

6  3  (l.OOOOOOOOOOOOOOO.OOOOOOOOOOOOOOOOE-KXX)) 

1  4  (3.589486770491939E-003,O.OOOOOOOOOOOOOOOE-K)00) 

2  4  (O.OOOOOOOOOOOOOOOE-KIOO.O.OOOOOOOOOOOOOOOE-KXX)) 

3  4  (0.000000000000000E-K)00,0.000000000000000E-K)00) 

4  4  (0.000000000000000E-K)00,0.000000000000000E-K)00) 

5  4  (-2.334981756O2S643E-OO2,O.0O0O0OO0OOOO0OOE-K)OO) 

6  4  (O.OOOOOOOOOOOOOOOE-KXIO.O.OOOOOOOOOOOOOOOE-KXX)) 
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initial  eigoivector  matrix  (coot) 


row  column  ( real ,  imaginaiy ) 

1  5  (-7.539018574755219E-002,1.439135915971290E-012) 

2  5  (1.405590215265010E-011.1.021808214943804E-001) 

3  5  (-3.808102176268219E-017.-2.617711356936501E-018) 

4  5  (-1.757523295343237E-010.-6.646921145391658E-001) 

5  5  (l.OOOOOOOOOOOOOtXi.O.O^^ 

6  5  (-1.999332831  133622E-017,4.141361260090141E-017) 

1  6  (8.609626586666082E-017.1.470913273776291E-017) 

2  6  (-9.403731537416513E.017,4.1995e3329593159E-016) 

3  6  (3.092693812092022E-014,6.942300233683507E-002) 

4  6  (1.391287225318998E-015.-6.660459595692267E-015) 

5  6  (-1.332598592779162E-015.-3.703060854033629E-016) 

6  6  (l.OOOOOOtKXKWOOOOAOO^ 


The*/ matrix; 

row  column  ( real ,  imaginary ) 

1  1  (O.OOOOOOOOOOOOOOOE+OOO 

2  1  (O.OOOOCXKIOOOOOOOOE+OOO.O.OOOO^^ 

3  1  (0.0()(K)()()()0()()0()000EHK)()0,0.()^^ 

4  1  (O.OOOOOOOOOOOOOOOE-KWO.O.OOOC^ 

5  1  (O.OOOOOtXXlOOOOOOOE+OOO.O.OOOO^^ 

6  1  vOOOOOOOOOOOOOOOOE+OOO.O.OOOOOtK^^ 

1  2  (0.0()0()()()()()00()OOOOE4<)00,0.0000^ 

2  2  (0.0(K)()0()()()(K)()()000E-K)003.853941 8253070^ 

3  2  (O.OOOOOCKXKKKKKIOOE+OOO.O.OOOOOO^^ 

4  2  (O.OOOOOOOOOOOOOOOE+OOO.O.OOOOOO^^ 

5  2  (0.(X)()()()()()0(X)00()00E4<)00,0.000()0(^^ 

6  2  (0.00()(K)()()00000(X)0E4<)()0,0.0000()0()(^^ 

1  3  (0.()(X)()00()()()()0()000E4<)()0,0. 

2  3  (0.0()()00()()(X)()()()()()0E4<)(X),0.(W 

3  3  (O.OOOOOCXXWOOOOW 

4  3  (0.0()()()0()0()000()0()0E4<)00,0.00(X^ 

5  3  (0.()()0C)00()0<)00()()()0E4<)00,0.()(X)^^ 

6  3  (O.OOOOOOOOOOOOOOOE+OOO.O.OOOCXlOOt^ 
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J  matrix  (cont.) 


row  column  ( real ,  imaginary ) 

1  4  (l.OOOOOOOOOOOOOIXi.O.OOOOOOOOlM 

2  4  (0.()0()0()()()0()0(X)0()OE+<)00.0.()000()0^^ 

3  4  (0.000()0()()0()()()0()00E-K)00,0.()000(X)(X)()^^ 

4  4  (O.OOOOOOOOOOOOOOOE^OOO.O.OOOO^ 

5  4  (0.()<)()()(X)()0()00()(XX)E4<)00.0.()00()()0^^ 

6  4  (O.OOOOOOOCKXiOOOOOE+OOO.OW^ 

1  5  (0.()()()0()()(X)(X)()()()()0E4<)00,0.0a 

2  5  (O.OCXWOOOOOOOOOOOE+KiOO.O.OOOOOOO^ 

3  S  (0.()()(X)()()()00()0()()()OE4<)00.0.0()0()0^ 

4  5  (O.OOOlXWOOOOOOOOOE+iXiO.O.OCKKK^^ 

5  5  (0.()()()()0()()00()()()()00E+<)()0,-8.85394182 

6  S  (0.0()()()0()()()<)00()(K)0E4<)()0.0.00()(X^^ 

1  6  (O.OOOOOOOOOOOOOOOE+OM.O.OOW 

2  6  (0.()0()()()()()()000()(K)OE4<)00,0.000()(^^ 

3  6  (O.lXlOOOOOOOOOOOOOE^KiOO.O.lX^ 

4  6  (0.()()()00()000<KX)()()0E4<)()0,0.()0()(X)^^ 

5  6  (0.()()()()()()()()00()()()<)OE4<)00,0 

6  6  (0.0()()()()()()()()0()()()()OE+<^ 
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Appendix  B:  Ephemerides  from  February  1967 
The  Mowing  data  was  collected  from  the  i^mencon  Bphemeris  ondNmtiical 
Almanac  (U  S.  Naval  Obsavatofy,  1%7:19,53)  for  Februaiy  1%7: 


Date/Time  Moon's  Apparent  Longitude  Moon’s  Apparent  Latitude 

O  I  W  O  I  N 

7 Feb/ 0000 UT  289.d0  IS.dO  55.35d0  -4.d0  38.d0  35.46d0 

8 Feb/ 0000 UT  301.d0  56.d0  55.63d0  -4.d0  56.d0  09.13d0 

9  Feb/0000  UT  314.d0  27.d0  13.68d0  -4.d0  59.d0  07.25d0 

10 Feb/ 0000 UT  326.d0  46.d0  52.49d0  -4.d0  47.d0  57.13d0 

1 1  Feb/ 0000 UT  338.dO  56.d0  15.14d0  -4.d0  23.d0  42.42d0 


Date/Time  Moon's  Horizontal  Parallax  Sun's  Longitude 

O  «  M  O  I  N 

7  Feb/0000  UT  O.dO  56.d0  01.692d0  317.d0  31.d0  10.3d0 

8 Feb/ 0000 UT  O.dO  55.dO  34.266d0  318.d0  31.d0  58.3dO 

9 Feb/ 0000 UT  O.dO  55.dO  09.022d0  319.d0  32.d0  45.2d0 

10 Feb/ 0000 UT  O.dO  54.d0  46.509d0  320.d0  33.dO  30.8d0 

1 1  Feb/ 0000 UT  O.dO  54.d0  27.580d0  321.d0  34.d0  IS.OdO 
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Appendix  C:  Ephemerides  from  February  1986 
This  appendix  contains  the  information  that  was  used  to  run  the  alternate  new  moon 
data  for  verification  of  the  errors  in  the  preliminary  results.  The  following  data  was  collected 
from  The  Astronomical  Almanac  (U.S.  Naval  Observatory,  1986;C4,D6)  for  February  1986: 


DateA'ime 

Moon's  Apparent  Lxxigitude 

0 

Moon's  Apparent  Latitude 

O 

7Feb/0000UT 

291.38d0 

-4.91d0 

8Feb/0000UT 

305.52d0 

-5.01d0 

9Feb/0000UT 

319.44d0 

-4.81d0 

10  Feb /OOOOUT 

333.07d0 

-4.34d0 

11  Feb/0000  UT 

346.35d0 

-3.64d0 

Date/Time 

Moon's 

O 

Horizontal  Parallax 

«  H 

Sun's 

O 

Longitude 

1  N 

7Feb/0000UT 

O.dO 

59.d0 

05.18d0 

317.d0 

55.dO 

01.47d0 

8  Feb /OOOOUT 

O.dO 

58.dO 

41.18d0 

318.d0 

55.dO 

49.69d0 

9  Feb /OOOOUT 

O.dO 

58.d0 

09.23d0 

319.d0 

56.d0 

36.72d0 

10  Feb /OOOOUT 

O.dO 

57.d0 

31.30d0 

320.d0 

57.dO 

22.47d0 

11  Feb/0000  UT 

O.dO 

56.d0 

50.16d0 

321.d0 

58.d0 

06.84d0 

The  new  moon  selected  from  1986  occurred  9  February  1986  at  0055  UT,  tiierefore  the 
interpolation  variables  and  subsequent  initial  state  vector  for  tiiis  case  are: 


A  =1.0  (dtQrs)  a  =  (days) 

1440.0 

u_2  =  7  Feb  OOOOUT 
tt_,  =8Feb  OOOOUT 
Mo  =  9  Feb  OOOOUT 
M,  =10  Feb  OOOOUT 
Mj  =llFeb  OOOOUT 


Jc(0)  =  [-0.99748 


-6.9847E-7  -2,1084£--4 


1. 575  !£•  -  3  -0. 96495  7. 7049£:  -  4]^ 
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